/* 
 * File:   UpgmaAlgorithm.cpp
 * Author: loskana
 * 
 * Created on March 5, 2011, 2:33 PM
 */

#include <cstdlib>
#include <iostream>
#include <math.h>
#include <stdio.h>

#include "Align.h"
#include "ScoringS2S.h"
#include "SequenceData.h"
#include "NWAlign.h"
#include "AGPFunction.h"

#include "UpgmaAlgorithm.h"
#include "NewickTreeFactory.h"

namespace Biopool
{

#if !defined(max)
#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) < (b) ? (a) : (b))
#endif

UpgmaAlgorithm::UpgmaAlgorithm(std::vector<std::string> & inputSequences) :
sequences(inputSequences) {
	
	// set distances parameters
	this->distanceGapPenalty = 10;
	this->distanceGapExtension = 2;
	
	// create Align classes 
    std::ifstream matrixFile("data/blosum62.dat");
	if (!matrixFile) 
        throw "Error opening structure matrix file.";
        
    this->subMatrix = new SubMatrix(matrixFile);
    this->gapFunction = new AGPFunction(this->distanceGapPenalty, this->distanceGapExtension);
}

UpgmaAlgorithm::UpgmaAlgorithm(const UpgmaAlgorithm& orig) :
sequences(orig.sequences) {
}

UpgmaAlgorithm::~UpgmaAlgorithm() {
	delete(this->subMatrix);
	delete(this->gapFunction);
}

// TODO: fa ritornare string con serializzazione
NewickTree UpgmaAlgorithm::execute() {

    std::cout << "\t \t \t * Starting Upgma Algorithm... * " << std::endl << std::endl;

    // Inizialization
    
    std::cout << "\t \t \t * Calcolo delle distanze iniziali... * " << std::endl << std::endl;

    // calculate distance matrix
    this->calculateInitialDistanceMatrix();

    // create initial clusters
    for (unsigned int i = 0; i < sequences.size(); i++) {
        intSet * setTemp = new intSet;
        setTemp->insert(i);
        NewickTree::Node * temp = new NewickTree::Node(setTemp);
        this->clusters.push_back(temp);
    }

    int numClusters = this->clusters.size();
    NewickTree::Node * nodeRoot = 0;

    std::cout << "\t \t \t * Unione dei clusters... * " << std::endl << std::endl;
    
    while (numClusters > 1) {

        // a, b: index dei cluster da accoppiare in questa iterazione
        // l: distanza tra i cluster
        int a, b;
        float l;
        // find minimum L(a,b)
        this->getMinimumDistance(a, b, l);

        std::cout << "Unione: a" << a << ", b: " << b << ", ed l: " << l << "\n" << std::endl;

        this->updateDistance(a, b);

        nodeRoot = this->updateClusters(a, b, l);

        // aggiorna contatore dei clusters
        numClusters--;
    }
    
    std::cout << "\t \t \t * Costruzione dell'albero completato * " << std::endl << std::endl;

    NewickTree newickTree(nodeRoot);
    NewickTreeFactory treeFactory;
    return newickTree;
    
    // TODO: fa ritornare string con serializzazione
    std::cout << "\t \t \t * Serializzazione del newick tree...  * " << std::endl << std::endl;
        
    // return treeFactory.serialize(newickTree);

}

float UpgmaAlgorithm::calculateDistance(std::string firstSeq, std::string secondSeq) {

    this->sequenceData = new SequenceData(2, firstSeq, secondSeq);
    this->scoringScheme = new ScoringS2S(this->subMatrix, this->sequenceData, 0, 1.00);
    this->align = new NWAlign(sequenceData, this->gapFunction, this->scoringScheme);

	this->align->doMatchScore(cout);
	
    // get alignment
    vector<string> match = this->align->getMatch(); 
    cout << std::endl << "Match: " << std::endl;
	for (unsigned int i = 0; i < match.size(); ++i)
        cout << match[i] << std::endl;
        
    firstSeq = match[0];
    secondSeq = match[1];
    
    // compute residues that differs
    int diffResidues = 0;
    int n = min(firstSeq.size(), secondSeq.size());

    for (int i = 0; i < n; i++) {
        if (firstSeq[i] != secondSeq[i])
            diffResidues++;
    }

    float residuesFract = ((float) diffResidues) / n;

    if (residuesFract > 0.75)
        residuesFract = 0.75;
        
    delete(this->sequenceData);
    delete(this->scoringScheme);
    delete(this->align);

    return residuesFract;
}

int UpgmaAlgorithm::getDistanceIndex(int firstId, int secondId) {

    int n = this->sequences.size();

    if (firstId > n - 1 || firstId < 0)
        throw "Accesso alla distanceMatrix con firstId errato";

    if (secondId > n - 1 || secondId < 0)
        throw "Accesso alla distanceMatrix con firstId errato";

    if (firstId == secondId)
        return 0;

    int i = max(firstId, secondId);
    int j = min(firstId, secondId);

    return (j + ((i - 1) * i) / 2);
}

void UpgmaAlgorithm::calculateInitialDistanceMatrix() {

    std::cout << "Calculating Initial Distance Matrix... " << std::endl;

    int n = this->sequences.size();
    // dimensione della matrice delle distanze, ottimizzata per l'utilizzo della sola upper matrix senza diagonale
    int m = ((n * n) - n) / 2;

    this->distanceMatrix.resize(m);

    // calcola la matrice di distanza - solo upper matrix, senza diagonale)
    for (int i = 0; i < n; i++) {
        // DEBUG std::cout << "i: " << i << "\n" << std::endl;

        for (int j = 0; j < i; j++) {

            // DEBUG std::cout << "j: " << j << "\n" << std::endl;

            // calculate distance between sequence i and j
            float residuesFract = this->calculateDistance(
                    this->sequences[i], this->sequences[j]);
            float distance = -3.0 / 4 * log(1 - 4.0 * residuesFract / 3);

            printf("   -> distance %d e %d : %4.2f \n ", i, j, distance);

            // inserisco il valore nella matrice, con allocazione efficiente
            int index = this->getDistanceIndex(i, j);
            this->distanceMatrix[index] = distance;
        }
    }
}


// cerca a, b , cluster con distanza minima

void UpgmaAlgorithm::getMinimumDistance(int & a, int & b, float & l) {

    // TODO?
    l = 100000;
    int n = this->sequences.size();

    for (int i = 0; i < n; i++)
        for (int j = 0; j < i; j++) {
            int index = this->getDistanceIndex(i, j);
            float distance = this->distanceMatrix[index];
            // se viene trovata una distanza migliore
            if (distance > 0 && distance < l) {
                a = i;
                b = j;
                l = distance;
            }
        }

    if (!l > 0) {
        throw "Non è stata trovata nessuna distanza maggiore di 0..";
    }
}

NewickTree::Node * UpgmaAlgorithm::updateClusters(int a, int b, float l) {

    // prendi i cluster a e b, crea un nuovo nodo che li contiene
    NewickTree::Node * nodeA = this->clusters[a];
    NewickTree::Node * nodeB = this->clusters[b];

    // aggiungi al set di K tutte le sequenze dei cluster di a e b
    intSet * seqK = new intSet;
    intSet::iterator setIt;
    intSet * seqA = nodeA->getSequencesSet();
    for (setIt = seqA->begin(); setIt != seqA->end(); setIt++)
        seqK->insert(*setIt);
    intSet * seqB = nodeB->getSequencesSet();
    for (setIt = seqB->begin(); setIt != seqB->end(); setIt++)
        seqK->insert(*setIt);

    // crea  il nuovo cluster K, ad altezza l/2, con setK, e figli A e B
    NewickTree::Node * nodeK = new NewickTree::Node(l / 2, seqK, nodeA, nodeB);

    // for (setIt = seqK.begin(); setIt != seqK.end(); setIt++)
        // DEBUG std::cout << "--- >" << *setIt << "\n" << std::endl;

    // aggiorna vettore dei nodi
    this->clusters[a] = nodeK;
    this->clusters[b] = 0;

    return nodeK;
}

// aggiorna i costi della tabella distanze
// in a, mette i conti di tutte le sequenze verso (a,b)
// in b, mette tutto a -1

void UpgmaAlgorithm::updateDistance(int a, int b) {

    // DEBUG std::cout << "updateDistance" << "\n" << std::endl;

    int sizeA = this->clusters[a]->getSequencesSet()->size();
    int sizeB = this->clusters[b]->getSequencesSet()->size();

    int n = this->sequences.size();

    // aggiorna costi: nuovo cluster A
    for (int j = 0; j < a; j++) {

        int indexAL = this->getDistanceIndex(a, j);
        int indexBL = this->getDistanceIndex(b, j);
        this->distanceMatrix[indexAL] = 
                ((this->distanceMatrix[indexAL] * sizeA +
                this->distanceMatrix[indexBL] * sizeB) / (sizeA + sizeB));
        // DEBUG std::cout << "metti in " << indexAL << "valore" << this->distanceMatrix[indexAL] << "\n" << std::endl;
    }

    for (int i = a; i < n; i++) {

        int indexAL = this->getDistanceIndex(a, i);
        int indexBL = this->getDistanceIndex(b, i);
        this->distanceMatrix[indexAL] =
                ((this->distanceMatrix[indexAL] * sizeA +
                this->distanceMatrix[indexBL] * sizeB) / (sizeA + sizeB));

        // DEBUG std::cout << "metti in " << indexAL << "valore" << this->distanceMatrix[indexAL] << "\n" << std::endl;
    }

    // aggiorna costi: vecchio cluster B
    for (int j = 0; j < b; j++) {
        int index = this->getDistanceIndex(b, j);
        this->distanceMatrix[index] = -1;
    }

    for (int i = b; i < n; i++) {
        int index = this->getDistanceIndex(b, i);
        this->distanceMatrix[index] = -1;
    }

}

} /* namespace */
